Projected wave functions for fractionalized phases of quantum spin systems 
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Gutzwiller projection allows a construction of an assortment of variational wave functions for 
strongly correlated systems. For quantum spin S = 1/2 models, Gutzwiller-projected wave functions 
have resonating-valence-bond structure and may represent states with fractional quantum numbers 
for the excitations. Using insights obtained from field-theoretical descriptions of fractionalization in 
two dimensions, we construct candidate wave functions of fractionalized states by projecting specific 
superconducting states. We explicitly demonstrate the presence of topological order in these states. 



I. INTRODUCTION 

In the last several years, considerable theoretical 
progress has been made in understanding the possibil- 
ity of fractional quantum numbers for the excitations of 
strongly interacting systems. A particular focus of atten- 
tion has been on realizing such phases in quantum spin 
systems in two spatial dimensions pi. In conventional 
quantum phases of spin- 1/2 magnets, such as an anti- 
ferromagnet or a spin-Peierls phase, the spin excitations 
carry quantum numbers S = 1 or higher multiples. In 
contrast, in a fractionalized phase, there are excitations 
(dubbed spinons) that carry spin S = 1/2. Various kinds 
of evidence for the presence of such fractionalized phases 
in spin models on a number of diverse lattices have been 
presented in Ref. 0. A detailed theoretical study |^-||] 
of fractionalization in two dimensions has revealed that 
in the simplest cases such spinon excitations are neces- 
sarily accompanied by gapped "vortex" -like excitations: 
visons. The visons carry no spin but have a long-range 
statistical interaction with the spinon: When a spinon is 
taken all the way around a vison, the wave function of 
the system changes sign. Thus the visons have an Ising- 
like character: two visons can annihilate each other and 
are equivalent to no vison at all. A concise mathematical 
description of this long ranged statistical interaction is 
given by assigning a Zi gauge charge to the spinons and 
a corresponding Zi gauge flux to the visons. 

A precise theoretical characterization |||| of fraction- 
alized phases serving to distinguish them from more con- 
ventional phases is provided by the notion of topological 
order — a concept elucidated by Wen in his work 
on the quantum Hall effect. When placed on a manifold 
with a non-trivial topology, a fractionalized phase has a 
number of locally similar but globally distinct states that 
differ by whether or not a vison threads each hole of the 
manifold. 

If a spin model possesses a fractionalized phase, what 
does its wave function look like in terms of the origi- 
nal spins? In this paper, we argue that the wave func- 
tions describing a fractionalized state may be obtained 
from the original RVB construction of Anderson pj . We 
consider states obtained by projecting the wave function 



of a BCS superconductor (on a lattice) to the Hilbcrt 
space of exactly one electron per site (this procedure 
is also known as Gutzwiller projection [0). Differ- 
ent topological sectors are obtained by projecting states 
with/without superconducting vortices (or, equivalently, 
with different boundary conditions). We motivate this 
construction from the field-theoretic description of frac- 
tionalized phases. 

Projected wave functions have long been studied 
|l2| |4j to gain insight into the properties of the cuprate 
materials. However, as we show below, it is not guaran- 
teed that the resulting wave function describes any frac- 
tionalized phase. The state that has been the subject 
of most attention is the "d-wave RVB state" obtained 
by projecting a BCS d x i^ y i superconductor with near- 
est neighbor hopping and nearest neighbor pairing in- 
teractions. We argue below that such a superconduct- 
ing state has some special symmetry properties due to 
which it is not expected to lead to a fractionalized state 
after projection (at half- filling) . We confirm this by a 
direct numerical calculation on this state. It is not clear 
whether the nearest neighbor d-wave RVB state repre- 
sents a good trial wave function for any stable phase 
of a spin system. Addition of a number of perturba- 
tions, such as for instance next-nearest-neighbor hopping, 
to the nearest-neighbor ei-wave superconductor removes 
the special symmetry properties — the resulting state 
is then expected to lead, after projection, to a fraction- 
alized state. We demonstrate explicitly the topological 
order in such projected states. A closely related state 
was recently considered in Ref. jl5| at finite doping and 
was argued to provide a good description of cuprate phe- 
nomenology. We note that recent experiments []l6 17 
place constraints on the applicability of fractionalization 
ideas to the cuprates. Theoretical description of the 
cuprates with wave functions corresponding to fraction- 
alized states must hence be performed with caution and 
examined for consistency with experiments. 

Recently, a particular projected state has been argued 
to represent the ground state of the J\ — Ji Heisenberg 
model on the square lattice p8| . We argue that such a 
state is also expected to have topological order, and hence 
be fractionalized. We explicitly construct four globally 
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distinct states on the torus to demonstrate the topologi- 
cal order. 

If the projection does indeed lead to a fractionalizcd 
state, then we may construct wave functions for its exci- 
tations as follows. The distinct excitations of a supercon- 
ductor are the BCS quasiparticles and the vortices. As 
originally suggested by Anderson jjj, projection of the 
BCS quasiparticles directly leads to the spinons. We ar- 
gue that the visons may be obtained by projecting the 
wave function of the BCS hc/2e vortex. This observation 
may be directly exploited to check for topological order in 
the projected wave functions. Consider putting the sys- 
tem on a torus. We may project the four superconduct- 
ing states obtained by threading or not threading a hc/2e 
vortex through each hole of the torus. If the projected 
state is fractionalized, and hence topologically ordered, 
these four states must be (after projection) orthogonal to 
each other in the thermodynamic limit though they are 
locally similar. 



II. GENERALITIES ON GUTZWILLER 
PROJECTION 

The original point of view on the Gutzwiller projected 
wave function || is that it describes a state obtained by 
quantum disordering a Heisenberg Neel antiferromagnet. 
However, as Gutzwiller projection (at half filling) freezes 
out the charge fluctuations inherent in the unprojected 
superconducting state, it may equally well be viewed as 
describing a state obtained by quantum disordering a su- 
perconductor. Recent theoretical work (|J^] has sought 
to understand the properties of Mott insulators by view- 
ing them as quantum disordered superconductors. In this 
approach, insulating phases are regarded as condensates 
of the vortices of the superconductor. A fractionalized 
Mott insulator is obtained by pairing and condensing the 
BCS vortices while the unpaired hc/2e vortex remains 
gapped. At energy scales well below the charge gap, the 
excitations in such a Mott insulator are the spinon and 
the vison. The spinons are precisely what become of the 
fermionic BCS quasiparticles once the hc/e vortices con- 
dense. The vison on the other hand is the remnant of the 
unpaired hc/2e vortex - due to the hc/e condensate, this 
survives with only a Z2 character. 

It is clear that this point of view fits in nicely with 
the idea of Gutzwiller projecting superconducting states 
to obtain wave functions for Mott insulators. However 
not all choices of the unprojected superconducting state 
are guaranteed to lead to fractionalized states on projec- 
tion (see Sections below). In this section, we discuss the 
nature of projected vortex states and their use in con- 
structing the wave function of the vison. 



A. RVB wave functions 

We begin by quickly reviewing the interpretation of the 
Gutzwiller-projected superconducting state as a resonat- 
ing valence bond wave function. Consider the ground 
state of a BCS Hamiltonian on a L x L square lattice (L 
even) : 

ij a 

Its wave function has the form 

\BCS)=Y[(u- k + v j :cy_ h )\0) (2) 

k 

in standard notation. We restrict attention to spin- 
singlet, time-reversal-invariant ground states that also 
preserve all the symmetries of the lattice. The vectors 
k take values in the first Brillouin zone. We assume 
periodic boundary conditions in both spatial directions 
so that k = 2 ^-(m x ,m y ) with m x ,m y integers (the sys- 
tem may be thought of as residing on a torus). We may 
rewrite this wave function as 

\BCS) oc exp f £ 0®4t C U) |0) (3) 

with g(k) = Vz/ur. So long as g(k) = g(—k), this state 
is guaranteed to be a spin singlet. 

It will be useful to consider a first quantized version 
^({xi, Xi}) obtained by fixing the total number of parti- 
cles to be L 2 . The (a?j, x%) are the position and spin state 
of the i'th electron. We have 

The wave function ^({xi, X i}) is single- valued on the 
torus. 

The Gutzwiller-projected state \RVB) is obtained by 
projecting W further onto a subspace where there is ex- 
actly one particle per lattice site. The resulting wave 
function may be written as a sum over various valcncc- 
bond configurations on the lattice: 

\RVB) = P G \BCS) oc J2 A "bc \vbc) (5) 

vbc 

where \vbc) denotes any particular valence bond covering. 
A v bc is the corresponding amplitude and is given by the 
product of g(xi — Xj) over all valence bonds connecting 
sites (ij) appearing in the state \vbc). Here g(xi — xj) 
is the Fourier transform of the function g(k) introduced 
above. 
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B. BCS vortices and their projections 

If the ground state of a spin system is correctly de- 
scribed by Gutzwiller projecting a BCS state, it is natural 
to construct excitations of the spin system by similarly 
projecting the excitations of the BCS superconductor. 
As argued by Anderson, projecting the BCS quasipar- 
ticle state leads to a wave function for neutral spin- 1/2 
spinons. The other distinct excitations of the supercon- 
ductor are the vortices. What happens to these under 
the Gutzwiller projection? To answer this question, con- 
sider a wave function of a superconductor describing a 
state where a nhc/2e vortex threads the torus along the 
x-direction. If n is even, the (first-quantized) wave func- 
tion may be taken simply as 



(6) 



It is readily verified that this is a spin-singlet wave func- 
tion. The first quantized hc/2e vortex wave function 
(with L 2 particles) may now be obtained straightfor- 
wardly as 



*' « m,xi}\ £</(f,ao4(*)4^) 



with 



g\x,x?) = e i ^ + ^g a . p {x-x!), 



L 2 /2 



10} (10) 



(11) 



= cos (#.(af- a*)) (12) 



where 9i = is the angular coordinate of the i'th par- 
ticle. Thus, the vortex wave function for even n may be 
simply obtained by multiplying the ground state wave 
function by a phase factor (which depends on the con- 
figuration of the particles). Upon Gutzwiller projection, 
there is one particle at every site; consequently the pro- 
jected vortex state (for n even) is trivially related to the 
projected ground state by an overall phase factor. Thus 
the projected even vortex state does not lead to a distinct 
state of the spin system. 

For n odd, the vortex wave function in the supercon- 
ductor is not obtained by simply multiplying the ground 
state by a phase factor. Consider n = 1. The naive guess 



(7) 



violates the physical requirement that all legitimate elec- 
tron wave functions must be single-valued on the torus. 
The correct wave function is constructed by multiply- 
ing by the phase factor above, and replacing f by a 
wave function that satisfies antiperiodic boundary con- 
ditions along the y-direction (and periodic along the x- 
dircction). Such a replacement may be obtained by con- 
sidering the ground state of H with antiperiodic bound- 
ary conditions along the y-direction: 



(8) 



with k' = k + j-y and k = ^{m x ,m y ) as before. Multi- 
plying by the phase factor above has the effect of shifting 
all the momenta by J- along the y-direction. A legitimate 
hc/2e vortex state may therefore be constructed as 

l^> = n(«&+« & c^ T cl Si )|0> (9) 



Here, as before, gik') = The function g a . p {x — x') 
satisfies antiperiodic boundary conditions along the y di- 
rection, and periodic along the x direction: 



9a. P {r + Ly) = -g a . P {r) 
g a . P (r + Lx) = g a . p {r) 



(13) 
(14) 



This vortex state may now be projected to obtain the 
new RVB wave function given by 



\RVB') 



p g \bcs')^Ta: 



vbc 



vbc 



\vbc) 



(15) 



The amplitude A' vbc is given by the product of g'(xi-Xj) 
over all valence bonds connecting sites i and j appear- 
ing in the state \vbc). This wave function can be further 
simplified by noting that the factor exp(i|J J2 i y{) arises 
in every term on the right-hand side and is just a triv- 
ial multiplication by an overall phase factor, and hence 
may be dropped. (For an L x L square lattice with L 
even, this phase factor is simply equal to one). This then 
amounts to Gutzwiller projecting the state \BCS^ ) ob- 
tained by changing the boundary conditions on the un- 
projected state. The resulting amplitude for a valence 
bond (ij) is given by g a . P {xi ~ xj). 

Thus the projection of a hc/2e vortex potentially sur- 
vives as a non-trivial state in the spin system. In view 
of the discussion at the beginning of this section, it is 
clear that if the ground state described by \RVB) is in- 
deed fractionalized, then \RVB') describes a state with 
a vison threaded through along the x-direction. We will 
exploit this observation below to check for topological 
order in the projected wave functions. 
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C. Short ranged RVB states and dimer models 



III. TOPOLOGICAL ORDER: GENERAL 
CONSIDERATIONS 



It is instructive to make a short digression, and special- 
ize to superconducting states which are fully gapped, and 
where the function g(r) has a very short range of order a 
few lattice spacings. In that case, the resulting RVB wave 
function describes a state with only short ranged valence 
bonds (and presumably a full spin gap) . For such a state, 
it is expected that the function g(k) is smooth in fc-space 
so that we may approximate g{k') m g(k). In real space, 
this amounts to 



g a .p(x-x) « g(x 



x")cos[-(y 



y 



(16) 



In the limit of large L, the cosine factor is one every- 
where except for valence bonds that connect sites to the 
left of y = L to sites to the right of y — L for which it 
is —1. Thus the difference between the projected ground 
state and the projected hc/2e state may be summarized 
as follows: Consider drawing a vertical line parallel to 
the x-direction that cuts all the links of the lattice be- 
tween y — L and y = 1. In the projected vortex state, 
the amplitude for any valence bond that cuts this line ac- 
quires a factor (—1) relative to the ground state. In other 
words, the amplitude for any valence bond configuration 
where an odd number of valence bonds intersect this line 
is (—1) relative to the ground state. We note that this 
construction of the projected vortex wave function re- 
duces to the state considered by Read and Chakraborty 
in their pioneering early work 19 on topological order 
in RVB wave functions. Unfortunately, as already recog- 
nized in that paper | fl9[ it is not clear that the particular 
RVB state discussed there describes any stable phase of 
a spin system. 

An approximate description of a short ranged RVB 
state is through a quantum dimer model |^of . It is well- 
known |2p| , pl]] that the quantum dimer model on a torus 
has four distinct topological sectors classified by whether 
the number of dimers intersecting particular lines similar 
to the one introduced above is even or odd. Fractionaliza- 
tion in the original spin model is signalled by topological 
order in the quantum dimer model when the four topo- 
logical sectors become degenerate in energy. In this case, 
these even/odd sectors clearly correspond to the symmet- 
ric/antisymmetric combination of the state with no vison 
and a threaded vison. This is completely consistent with 
our construction of the vison wave function by projecting 
the hc/2e vortex. 

Our primary interest in this paper will be long-ranged 
RVB states where the spinons are gapless. We will ex- 
plore the properties of wave functions for such states by 
numerical calculations using the construction above. 



Though the projected hc/2e vortex survives as a non- 
trivial state of the spin system, it is still not necessary 
that the state \RVB) has topological order. The presence 
or lack thereof of topological order may be established by 
examining the following two conditions: 

(i) (RVB'\RVB) = as the system size goes to infin- 
ity, and 

(ii) The expectation values of all local physical opera- 
tors be the same in both states \RVB) and \RVB') again 
in the thermodynamic limit. 

The first condition is closely related to the presence of 
a vison gap in the bulk. Indeed, the states \RVB) and 
\RVB') differ by one vison tunneled across the cylin- 
der (or torus), and therefore (RV B\RV B') may be in- 
terpreted as the amplitude of such a tunneling event. 
The second condition guarantees that the distinction be- 
tween the states is not in any local properties but rather 
in global ones. In this section, we will examine some 
simple arguments that motivate the choice of particular 
superconducting wave functions which will lead to frac- 
tionalized states on projection. 

Consider the first condition. It is easy to check that the 
overlap (BCS'\BCS) goes to zero very rapidly as L be- 
comes large. This is expected physically as the vorticity is 
a good quantum number for the unprojected BCS state. 
However, this does not guarantee (RVB'\RVB) = due 
to the projection. To get better insight into what is 
needed for condition (i) to be satisfied, we will employ the 
following useful characterization of the Gutzwiller projec- 
tion. In the unprojected Hilbert space, the states at each 
site are |0) , |t) , |J.) , |t|) in obvious notation. Projection 
keeps only the states |T) , I!) at each site. At each site i, 
introduce the physical spin operator Si and the pseudo- 
spin operator Tf. 



o 1 t- 



T+ - c f c f 



(17) 

(18) 

(19) 
(20) 



Here T± = T?± Tf . The f; satisfy SU(2) commutation 
relations just like the physical spin operators Si. Fur- 
thermore, all components of Tj commute with all compo- 
nents of Si. Clearly the states |f) , ||) are singlets under 
the pseudo-spin rotation generated by the T il while the 
states |0) , |f |) form a doublet under the pseudospin ro- 
tation. Thus the Gutzwiller projection is equivalent to 
projection onto the singlet sector (at each site) of the 
pseudospin SU (2) rotation. This has the general impli- 
cation 
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Pg I ^unproj) = PgU | ^unproj) (21) 

for an arbitrary SU(2) rotation 

[/ = e ! £. 9 "*- f - (22) 

where the parameters 9i may be chosen independently 
on each site. Thus a local pseudospin SU(2) rotation 
of the unprojected state does not change the state after 
projection. We will therefore call this a gauge rotation. 

In view of the above, a sufficient condition for the or- 
thogonality of \RVB) and \RVB') is simply 

(BCS\U\BCS'} = (23) 

for any arbitrary pseudospin SU(2) rotation (We 
note that this is not in general a necessary condition for 
the orthogonality). 

The discussion has so far been general. Eqn. (|2^) im- 
poses some conditions on the nature of the BCS state 
which could lead to a wave function for a fractionalized 
state after projection. However Eqn. (|2^) is still not in 
a form which is directly useful in providing guidance in 
writing down such states. To get a more useful form, we 
specialize to a particular class of unprojected states. In 
general, the state \BCS) will be a linear superposition of 
states with different total particle number. On the other 
hand, the projected state \RVB) has exactly one parti- 
cle per site, and hence an exact total of L 2 particles on 
a lattice with L 2 sites. Now, for a general BCS state, if 
we plot the probability distribution of having a total of n 
particles, it will have a sharp peak at some average value 
no and will die rapidly for \n — no\ large. If no is differ- 
ent from L 2 , then the Gutzwiller projection picks out the 
tails of the original BCS wave function. In this case, the 
projected wave function may have very little to do with 
the unprojected one. For the projected state to retain the 
significant features of the spin physics of the unprojected 
state, it is clearly advantageous to require that the 
mean number of particles (before projection) is L 2 . If 
this is satisfied, the projection will have a gentler effect 
than if the mean number is different from one per site. 
To further soften the effect of the projection, it is clearly 
advantageous to require that the mean number be one 
for every unprojected state obtained from \BCS) by a 
pseudospin SU (2) gauge rotation. This ensures that the 
projection does not pick up the tails of the wave function 
in any gauge. 

The requirement that the average electron number on 
each site is one may be simply written as 

(If ) = (24) 

on all sites i. Further the requirement that this be true 
for every gauge-rotated ground state is equivalent to re- 
quiring that 

(Ti) = 0, (25) 



i.e. the average values of the generators of the gauge ro- 
tations vanish. 

From now on, we will specialize to unprojected states 
where the requirement Eqn. (|25| ) is satisfied. In this case, 
we may hope that we can use our intuition about the un- 
projected state to infer the properties of the projected 
state. When do we expect that Eqn. (|2^) will be sat- 
isfied for such a BCS state? Note that \BCS') differs 
from \BCS) in having a hc/2e vortex threaded through 
the hole of the cylinder. If the ground state of a physical 
system is given by the superconducting state \BCS), we 
can label it's excitations by their total vorticity quan- 
tum number. The state \BCS') has total vorticity ±1 
compared to \BCS), and hence is orthogonal. However, 
so long as U\ BCS) is also a superconducting state, it 
will also have a fixed vorticity which differs from that of 
\BCS') by an odd number. Consequently, in this case 
Eqn. (|2^) will be satisfied. If however for some U, the 
gauge rotated state U \BCS) is not a superconducting 
state, then it may be regarded as a coherent superpo- 
sition of electron wave functions each of which carry a 
definite vorticity. (The vorticity is not a good quantum 
number in a non-superconducting state). It is then pos- 
sible that \BCS') is not orthogonal to that particular 
gauge rotated state, so that Eqn. (|23| ) is not satisfied. 
At the very least, we lose the general reason requiring 
orthogonality. 

Based on the reasoning above, we make the following 
conjecture: 

An unprojected state where the constraint ( |25| ) is 
satisfied, and which cannot be rotated to a non- 
superconducting state by an SU{2) gauge rotation, will 
lead to a topologically ordered state after projection. 

This conjecture provides strong motivation for choos- 
ing particular superconducting states that we may 
project to obtain topologically ordered states. On the 
other hand, we do not provide any formal proof of this 
conjecture in this paper. It is supported by the reasoning 
in this section, and by our numerical results. 

Here we need to make a clear distinction between 
"superconducting" and "non-superconducting" states. 
Formally, we call a state "gauge-equivalent to non- 
superconducting" if it is invariant under a global 17(1) 
subgroup of the full group of gauge rotations (|2^). In- 
deed, a non-superconducting state is characterized by a 
fixed number of electrons. Then the U(l) group of global 
rotations by T z (i.e., the subgroup generated by the sets 
9i = 8z) produces only trivial multiplications by phase 
factors leaving the state invariant. We call such a state a 
"Z7(l) state" [or even a "SU(2) state" in the case when 
the maximal subgroup leaving the unprojected state in- 
variant is SU(2)]. In this terminology, a state that is su- 
perconducting in any gauge is called a " Z-i state" . This 



5 



classification is a sub-classification of a more detailed one 
introduced by Wen |2^| . 

A set of conditions for the U(l) invariance (being 
gauge-equivalent to non-superconducting) of a given 
state may be formally written in terms of projector op- 
erators onto occupied quasiparticle states. Consider the 
set of Bogoliubov-deGennes doublets (it£, Ug) participat- 
ing in the ground state (||). We can define a projector 
onto the linear subspace spanned by those doublets. In 
real space, this projector may be thought of as a set of 
2x2 matrices Py labeled by the site indices i and j and 
defined as 



p, = E 



uUj) v*(j)) 



(26) 



Under the SU(2) gauge transformation, these matrices 
transform as 



p. 



UiPijU] 



(27) 



For a non-superconducting state, all matrices Py- are si- 
multaneously diagonal. To see when a given state can 
be gauge rotated to such a non-superconducting state, it 
is useful to consider the chain products of such matrices 
starting and ending at the same point i: 



AiC Hr (PyP ifc ...P M ) 



(28) 



where (i, j, k, . . . , I, i) define a closed curve on the lattice. 
Ai [C] is a 2 x 2 matrix for each lattice site i and for each 
closed curve C starting and ending at that site. Now one 
easily verifies that all matrices Py may be simultaneously 
diagonalized if and only if for any i and for any pair of 
closed contours C and C both starting and ending at i, 



[AiP\,MC]]=0 



(29) 



Thus to check that a given wave function describes 
a Z2 state (i.e. cannot be gauge rotated to a non- 
superconducting form) , it is sufficient to verify that Ai [C] 
and Aj[C] do not commute for at least one choice of i, 
C, and C. 

Instead of checking whether the wave function may be 
rotated to a non-superconducting form, one may perform 
a similar test for the Hamiltonian (Q). The Hamiltonian 
can be rotated to a non-superconducting one (i.e con- 
taining no pairing terms) if and only if 



[B i [C],B i [C']} = 0, 



(30) 



for all i and C, where 



Bi[C\ = ll<- (HijHju ...H H ), II,, ( '.; X > 

(31) 

This is a convenient sufficient condition for being a 
U(l) wave function, but not a necessary one: in certain 



cases, a superconducting Hamiltonian may have a non- 
superconducting ground state (with a definite particle 
number). 

It is instructive to consider some specific examples of 
superconducting states to see how these conditions work 
in practice. Consider for instance a nearest neighbor d- 
wave superconductor where Uj , A^- are non-zero only on 
nearest-neighbor bonds, and take the values 



A; 



t 

±A. 



(32) 
(33) 



The plus sign in the second equation is for horizontal 
bonds, and the minus sign for vertical bonds. It is read- 
ily seen that the condition in Eqn. (25) is satisfied by this 
state. However, as is also easily seen, all the Bi\C] com- 
mute for this state. This is of course consistent with the 
well-known fact that this state can be gauge rotated to a 
pure hopping state (known as the staggered flux state). 
We therefore expect that this state will not lead to a 
fractionalized state after projection. This is confirmed 
by our numerical calculations below. 

Now consider adding the next-near-neighbor diagonal 
hopping t' to this state. The resulting Hamiltonian still 
describes a d x 2_ y 2 superconductor. However, a sim- 
ple calculation shows that there are two non-commuting 
Bi[C] matrices so that this can no longer be rotated to 
pure hopping. Addition of t' changes the mean den- 
sity away from one particle per site. This can however 
be compensated by adding an on-site chemical potential 
term to the unprojected Hamiltonian. The state thus 
constructed is therefore a good candidate for projecting 
to get a wave function for a fractionalized state. Below, 
in Section IV we demonstrate this by a numerical calcu- 
lation. 

Surprisingly, the criteria obtained above for the super- 
conducting state to describe a topologically ordered state 
after projection may also be motivated from a completely 
different point of view [^|J23|. Consider any spin S =1/2 
model in two dimensions. As is well-known, it is possible 
to use a representation of the spins in terms of fcrmionic 
spin-1/2 operators. This representation is exact so long 
as the constraint that the fermion occupation is one at 
each site is imposed on the Hilbert space. A popular ap- 
proach is to treat the resulting fermion Hamiltonian in 
mean field theory. At the mean-field level, the excitations 
are neutral spin-1/2 fermions described by a Hamiltonian 
of the general form Eqn. ([!]). The exact constraint inher- 
ent in the fcrmionic representation is replaced precisely 
by Eqn. (^5j) . This mean field state is capable of correctly 
describing a possible physical phase of the spin model as 
long as it is stable to fluctuations. As discussed by Wen 
p],p3"t , the criteria for the stability of the mean field state 
to fluctuations are precisely that expressed in Eqn. d3C|). 
If the mean field state is indeed stable to fluctuations, 
then we expect that the candidate wave function for the 
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physical state described by it is given by the Gutzwiller 
projection of the mean field wave function to the physical 
Hilbert space. 



IV. NUMERICAL RESULTS 



We further verify our conjecture numerically by test- 



ing the conditions (i) and (ii) of Section III on the square 
lattice in the toroidal geometry for several examples of 
the projected BCS wave functions. We label the ground 
states in the four topological sectors as |++), |H — }, 

| — +), and | ), according to the boundary conditions 

imposed in x and y directions. We employ the variational 
Monte Carlo method described in detail in Ref. jl3| ap- 
plied to the square LxL tori. 



t.A 





NND Dl D2 DD 

FIG. 1. The mean-field states generating the four wave 
functions NND, Dl, D2, and DD. The terms of the Hamil- 
tonian (Q) are shown on one plaquette of the lattice. The 
solid lines denote hopping, the dashed lines denote pairing, 
and the circles around vertices denote the chemical potential. 

We consider four different types of projected wave 
functions: the nearest-neighbor d x 2_ y 2 state and its three 
modifications by including hopping or pairing along the 
plaquette diagonals (Fig. |l]). All these wave functions 
may be parameterized by translationally invariant Hamil- 
tonians (|l|), and we can conveniently describe them in 
terms of the Fourier transform of Hij (defined above in 
Eq. @): 



H(k) 



£(k) A(k) 
A*(k) -C(-k) 



(34) 



The nearest-neighbor d x 2_ y 2-wave BCS state (further 
denoted as NND) is defined by its kinetic and pairing 
amplitudes: 



f(k) = -2(cosfc 

3 k x — cos ky) . 



A(k) = A (cos 



+ COS Ky ) 



(35) 
(36) 



In the projected NND wave function, the nearest- 
neighbor antiferromagnetic correlations arc maximized at 
the intermediate value of Ao « 0.55 Jl2|,^4| (the optimal 
values of Ao reported in these two references differ by 
several percent; the precise value of A is not important 
for our qualitative results concerning topological order). 
We find that the NND state after projection has no topo- 
logical order. This agrees with our conjecture, since the 
NND state can be rotated to the pure-hopping staggered- 
flux state by a SU(2) gauge rotation p|]2§]. 






10 12 14 16 18 



FIG. 2. Overlaps between wave functions with different 
boundary conditions on tori L x L as a function of the sys- 
tem size L. All wave functions have Ao = 0.55. (a): over- 
laps (+ — | — h). Stars: NND wave function. Circles: Dl 
wave functions (dashed line: t' = 0.5, [i = 0.417; solid line: 
t' — 1.0, fj, — 0.696). Squares: D2 wave functions (dashed 
line: t' = 0.3, n = 0.509; solid line: t' = 0.5, fj, = 0.794). 
Triangles: DD wave functions (dashed line: D xy — 0.4; solid 
line: D xy = 1.0). (b): overlaps for the Dl wave function with 
t' = 1.0, fi = 0.696. (c): overlaps for the D2 wave function 
with t' — 0.5, fi — 0.794. (d): overlaps for the DD wave 
function with D xy — 1.0. In plots (b)-(d), solid line corre- 
sponds to (H — | — h), dashed line — to (+ — dotted line 

— to (+ — | ). For these three wave functions, the overlap 

(++| ) is found to be zero for all system sizes, within the 

numerical accuracy (about 0.003). 

In addition to the NND state, we consider its three 
modifications which are Zi states (cannot be gauge ro- 
tated to a non-superconducting form) . In the first modi- 
fication (we denote it by Dl), we add hopping along one 
of the plaquette diagonals which amounts to replacing 
©by 

£(k) = — 2[cosfc x + cosfcy + t' cos(k x + k y )] — p . (37) 



The chemical potential /i is added for adjusting the aver- 
age particle density in the unprojected wave function in 
order to satisfy the mean- field constraint (|2^). 

The second wave function (dubbed D2) is analogous to 
the previous one, but with hopping along both plaquette 
diagonals: 



£(k) = — 2 [cos + cos k y 

+t'(cos(k x + k y ) + cos(k x 

The third wave function is d x 2_ 
tion proposed by Capriotti et al 



k v ))]-(x. (38) 



- d xy wave func- 
as a variational 
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ansatz for the J\—Ji Heisenberg model. This wave func- 
tion (denoted further as DD) has the nearest-neighbor 
form ( |35| ) of £(k), together with A(k) involving both 
nearest-neighbor pairing and pairing along the plaque- 
tte diagonals: 

A(k) = Ao(cos k x — cos k y + 2D xy sin k x sin k y ) . (39) 

One easily verifies that the projected DD wave function 
has all the symmetries of the square lattice, even though 
the unprojected wave function does not. This wave func- 
tion obeys the mean-field constraint (^5|) without a chem- 
ical potential term. 



0.35 




0.25 1 1 1 ' 1 ' 1 0.0 1 1 1 ' 1 1 ' 1 1 

4 6 8 10 12 14 16 4 6 8 10 12 14 16 18 20 



FIG. 3. Nearest-neighbor spin correlations S a p for wave 
functions with different boundary conditions on tori L x L 
as a function of the system size L. All wave functions have 
A = 0.55. (a): Dl wave function with t' = 1.0, fi = 0.696. 
(b): D2 wave function with t' = 0.5, ^ = 0.794. (c): DD 
wave function with D xy = 1.0. In plots (a)-(c): solid sym- 
bols, solid lines correspond to S+_; solid symbols, dashed 

lines — to S |_; empty symbols, solid lines — to S++; empty 

symbols, dashed lines — to S (d): the spread of local cor- 
relations AS as a function of the system size L. Solid circles: 
Dl wave function, t' — 1.0, fi = 0.696. Empty circles: Dl 
wave function, t' = 1.2, /i = 0.744. Solid squares: D2 wave 
function, t' = 0.5, n — 0.794. Empty squares: D2 wave func- 
tion, t' = 0.6, /i = 0.909. Solid triangles: DD wave function, 
D xy — 1.0. Empty triangles: DD wave function, D xy — 1.2. 
In all plots, the error bars are smaller than the symbol size, 
except for the rightmost points in plot (d) where the error 
bars are of the order of the symbol size. 

The first test for the topological order is the over- 
lap of the two wave functions on a torus projected from 
the mean-field states with different boundary conditions. 
Note that when nodes in the mean-field excitation spec- 
trum fall on points of the momentum lattice, the wave 
function becomes ill-defined. For our choice of lattice 
on the tori, this happens for the NND state where the 
nodes are fixed at (±-|,±^) points. Therefore, not all 



four topological sectors are realized for the NND state 
(for our choice of the lattice placement). However, the 
sectors +— and — h are well-defined for any L x L torus, 
and we take the overlap between these two wave functions 
as the first check of the topological order. Alternatively, 
this overlap may be viewed as the overlap of the wave 
function with its reflection with respect to the x-y diag- 
onal. 

We have calculated this overlap with the use of the 
variational Monte Carlo procedure p3| . The results for 
the NND wave function and for its three modifications 
Dl, D2, and DD, are presented in Fig. ^a. In all cases, 
we take Ao = 0.55; for Dl wave functions we take t' = 0.5 
(// = 0.417) and t' = 1.0 (// = 0.696); for D2 wave 
functions we take t' = 0.3 (// = 0.509) and t' = 0.5 
(/i = 0.794); for DD wave functions we take D xy = 0.4 
and D xy = 1.0. We see that the overlap for the NND 
wave function saturates at large system size (which im- 
plies the absence of topological order), while for the three 
other wave functions it decreases to zero with increasing 
system size. This is the first signature of the topologi- 
cal order in those wave functions. As a side remark, we 
mention that if we force fi = in the D2 wave function, 
the projected wave function has no topological order (the 
overlap saturates at large L), but such a wave function 
violates the mean- field-density constraint (p5|), and we 
do not discuss it here. 

For the wave functions deep in the "topologically or- 
dered" phase, we can also verify the orthogonality of all 

the four topological sectors ++, H — , — h, and . The 

four corresponding overlaps are plotted for each of the 
wave functions Dl, D2, and DD in Fig. ||b-d. All the 
overlaps for these wave functions decrease to zero with 
increasing system size, which indicates orthogonality of 
the four sectors. 

The second test for the topologically ordered state is 
the equality of local expectation values between differ- 
ent topological sectors (condition (ii) of Section |l|). We 
verify this condition for the nearest-neighbor spin-spin 
correlations. In the four topological sectors on the torus, 
there are four different nearest-neighbor correlators: 

S+- = -<+-k 2 ^+xl+->, 

S-+ = -<+-|«+yl + ->» 

S++ = -{++\a?a]\++}, (40) 
S- = -{— \a?a*\— ). 

In a topologically ordered state, the difference between 
these four quantities must rapidly decay with increasing 
system size. We compute these correlators numerically 
for the same three wave functions used before in the over- 
lap computations. The results are presented in Fig. ^a-c 
demonstrating that the four correlations indeed converge 
to a single value in large systems. This supports the 
contention that the distinction between the states is in 
global properties and not in local ones. 
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To quantify the rate of this convergence, we consider 

the mean square deviation of the four quantities , 

S |_, S++, and S In Fig. ||d, we plot this mean square 

deviation multiplied by the number of lattice sites 



AS: 



Ea/3 



Eq/3 S a fj 



1/2 



(41) 



as a function of system size. The finite-size effects are 
very strong because of the nodes in the spectrum. To 
clarify the size dependence of AS", we take one more wave 
function in each of the three classes Dl, D2, DD, with 
slightly different variational parameters (in addition to 
the wave functions considered previously). The data in 
Fig. |^d indicate that AS* remains approximately indepen- 
dent of system size. This corresponds to the difference in 
local correlations S a p decaying as L~ 2 with system size. 

This slow convergence of correlations among different 
topological sectors can be explained already at the mean- 
field level from the nodal singularities in the spectrum. 
At the nodal points, the vector (u^,v^) has a strong sin- 
gularity: as the wave vector k goes around the nodal 
point, the vector (u, v) rotates by half turn in the plane. 
This singularity translates into the L~ 2 dependence of 
local correlations on the boundary conditions (changing 
the boundary conditions amounts to shifting the lattice 
of vectors k). From our numerical results we see that this 
rate of convergence is preserved by the projection. This 
may serve as an indication that in the Zi states the pro- 
jection preserves the nodal spinons, as suggested by Wen 
p3[ . This L~ 2 law applies to the general case of local cor- 
relations. The energy expectation value is a very special 
local correlation which has a weaker singularity at the 
nodal point. At the mean- field level, the energy splitting 
between different topological sectors can be found to be 
L -3 (per site). This translates into the splitting in total 
energy decaying as L . We expect that this asymptotic 
law also holds after the projection. To verify this numeri- 
cally, one needs to compute the expectation values of the 
actual "energy" , which is a correlation function optimized 
for a particular Hamiltonian by adjusting the variational 
parameters. In this paper we do not identify spin Hamil- 
tonians for which our wave functions are optimal, and 
thus are not able to verify our expectation for the energy 
splitting. Note that the L _1 convergence in energy also 
means convergence of the energy of an individual vison, 
i.e., deconfinement of visons which is necessary for the 
topological order. 

We have also tested all the above classes of wave func- 
tions for the valence-bond crystal ( "spin-Peierls" ) order- 
ing |27|] . We have computed the correlations of the z- 
component of the singlet order parameter {S z (i)S z (i + 
x)S z (j)S z (j + x)) in systems as large as 18x18. We find 
this correlation function rapidly decaying with increas- 
ing the distance \i — j\. The decay is slowest for the D2 
wave function for which it appears to be close to R~ 2 



(for the t' = 0.5, /i = 0.794 wave function) or faster (for 
the t' = 0.6, n = 0.909 wave function). For Dl and DD 
wave functions the decay of such correlations is much 
stronger than in the D2 wave function, which indicates 
the absence of the valence-bond ordering. Note that the 
D2 wave functions for the values of parameters consid- 
ered in this paper exhibit a relatively large correlation 
length for the overlaps between topological sectors (see 
Fig. ^c for the data on the t' = 0.5, ji = 0.794 wave func- 
tion). Therefore one may expect that correlation func- 
tions have different behavior at larger and smaller length 
scale, and that our computations in relatively small sys- 
tems are therefore not completely reliable for determining 
the correct long-distance behavior of the correlations. A 
more detailed analysis of the valence-bond crystal order- 
ing and of its interplay with the topological order is left 
for future study. 



V. CONCLUSION 

In this paper, we have formulated conditions for the 
presence of topological order in RVB systems and have 
verified them for several specific examples of Gutzwiller- 
projected wave functions. Our results suggest that ap- 
propriate Gutzwiller-projected wave functions may rep- 
resent ground states of fractionalized phases of spin sys- 
tems. 

This work is only the first step towards describing the 
topologically ordered RVB states. For a better under- 
standing of the properties of the topological order, a 
more extensive quantitative study is needed. It should 
include an analysis of correlation lengths involved in the 
conditions (i) and (ii) for the topological order (in par- 
ticular, on cylinders/tori with different aspect ratios). 
Variational wave functions may provide a useful tool for 
studying quantum phase transitions between states with 
and without topological order. An extremely interesting 
question in this respect is the possibility of the coexis- 
tence of the topological order simultaneously with the 
antiferromagnetic Neel order or with the valence-bond 
crystal order. Such a coexistence of a topological and 
a conventional ordering may be tested by projecting ap- 
propriate superconducting states with the spin or trans- 
lational symmetry broken before projection. 

Of course, the study of the variational wave functions 
have physical implications only when the Hamiltonians 
are identified for which those wave functions are good 
trial states. Our test for the topological order may pro- 
vide a guidance in the search for microscopic spin Hamil- 
tonians that exhibit fractionalized ground states. 
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